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1  Introduction 


Forward  scattering  has  been  an  active  field  of  research  in  science,  mathematics,  and  engi¬ 
neering  over  the  past  several  decades  (see  e.g.  [3],  [4]).  The  most  straightforward  method  for 
the  solution  of  a  forward  scattering  problem  is  to  discretize  the  underlying  PDEs  directly, 
replace  the  derivatives  with  finite  differences,  and  solve  numerically  the  resulting  system  of 
linear  algebraic  equations.  However,  discretization  of  differential  equations  leads  to  matri¬ 
ces  with  high  condition  numbers,  with  the  attendant  loss  of  accuracy,  deterioration  in  the 
performance  of  iterative  methods,  etc.  Another  approach  is  to  convert  the  underlying  PDEs 
into  integral  equations  of  the  second  kind  (normally  referred  to  as  the  Lippmann-Schwinger 
equation),  discretize  the  latter  via  appropriate  quadrature  formulae,  and  deal  numerically 
with  the  resulting  linear  systems.  This  paper  constructs  a  class  of  high-order  quadrature 
formulae  applicable  to  the  Lippmann-Schwinger  equation  in  two  and  three  dimensions. 

1.1  Statement  of  the  Problem 

The  forward  scattering  problem  is  the  problem  of  determining  the  scattered  field  given  the 
parameters  of  the  scattering  structure  and  the  incident  field.  In  this  section,  we  formulate 
the  two-dimensional  forward  scattering  problem  for  the  Helmholtz  equation,  and  derive  the 
corresponding  Lippmann-Schwinger  equation. 

The  forward  scattering  problem  we  investigate  arises  from  the  time  domain  wave  equation 

d 2 

— ^(rr,  t)  =  c2[x )  •  V2^(x,  t),  (1) 

where  ^(x,  t )  is  the  value  of  the  scalar  field  at  a  point  x  at  time  t,  and  c{x)  is  the  local  speed 
of  wave  propagation  at  a  point  x.  In  order  to  solve  (1),  we  assume  that 

^(x,t)=Mx)eikcot,  (2) 

where  k  is  a  complex  number  with  non-negative  imaginary  part,  and  Co  is  the  speed  of  wave 
propagation  outside  of  the  scattering  structure.  Substituting  (2)  into  (1),  we  obtain 

(V2  +  k2)ipk(x)  =  k2V(x)ipk(x),  (3) 

where 

V(x)  =  1  -  (-A)2.  (4) 

c(x) 

Equation  (3)  is  the  well-known  Helmholtz  equation,  and  the  operator  (V2  +  k2)  is  known 
as  the  Helmholtz  operator.  For  any  point  x  outside  the  scattering  object,  c(x)  =  c0;  therefore, 
V(x)  =0  outside  the  scattering  object.  We  represent  the  field  ^k{x)  at  a  point  x  as  a  sum 
of  two  parts:  the  incident  field  ^(x)  and  the  scattered  field  'ifj'lcat(x),  i.e., 

Mx)  =  WW  +  ^ka\x).  (5) 

The  incident  field  satisfies  the  homogenous  Helmholtz  equation 

(V2  +  k2)^in(x)  =  0,  (6) 
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in  some  open  region  in  R2  containing  the  scatterer;  the  scattered  held  satisfies  the  Sommerfcld 
radiation  condition 

lim  -iki>r‘(x))=0.  (7) 

|x|^oo  O  |X| 

Combining  equations  (3),  (5),  and  (6),  we  obtain  the  equation  for  the  scattered  held 

(v2  +  k2)  rr\x)  -  e  v(X)  rkcat(x)  =  e  vw  ww.  (s) 

In  this  paper,  we  view  the  equation  (8)  with  -0^ cat  satisfying  the  Sommerfcld  condition  (7) 
as  the  principal  formulation  of  the  forward  scattering  problem.  A  standard  approach  to  the 
numerical  solution  of  (8)  is  to  convert  (8)  into  the  well-known  Lippmann-Schwinger  equation, 
which  is  an  integral  equation  of  the  second  kind,  as  follows  (see,  for  example,  [5]). 
Convolving  (8)  with  a  Green’s  function  for  the  equation 

(V2  +  k2)Gk(x,y)  =S(x-y),  (9) 

we  obtain 

rkcat(x)-k2  [  Gk(x,y)V(y)rkcat(y)dy  =  k2  [  Gk(x,y)V(y)W(y)dy,  (10) 

J  D  J  D 

which  is  an  integral  equation  of  the  second  kind;  in  (10)  above,  D  denotes  the  region  in  space 
where  the  scatterer  is  located.  As  is  well-known,  in  two  dimensions,  the  Green’s  function 
Gk(x,y)  satisfying  the  condition  (7)  is 

Gk(x,y)  =  -l-H^\k\\x-y\\),  (11) 

where  H^\k\\x  —  y||)  is  the  Hankel  function  of  the  hrst  kind  of  order  zero.  We  will  define 
the  operator  L  :  L2(D)  —>  L2(D)  by  the  formula 

L(tp)(x)=  [  Gk(x,y)V(y)'ip(y)dy,  (12) 

J  D 

and  observe  that  a  large  part  of  this  paper  is  devoted  to  the  construction  of  accurate  dis¬ 
cretizations  of  L. 

1.2  Overview 

A  number  of  algorithms  exist  for  the  modeling  of  acoustic  scattering;  since  we  are  interested 
in  frequency  domain  results,  we  have  concentrated  on  frequency  domain  (as  opposed  to  time- 
domain)  models.  The  usual  approach  to  such  problems  is  to  convert  the  scattering  problem 
into  the  Lippmann-Schwinger  equation,  and  solve  the  latter  iteratively  (integral  equations  of 
the  second  kind  being  much  more  amenable  to  iterative  techniques  than  the  straightforward 
discretizations  of  underlying  partial  differential  equations  (PDEs)).  In  addition,  the  use  of 
the  Lippmann-Schwinger  equation  obviates  the  need  to  impose  the  radiation  (Sommerfeld) 
condition  at  the  boundary  of  the  grid,  since  the  “background”  Green’s  function  (11)  imposes 
the  Sommerfeld  condition  automatically. 
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Historically,  there  have  been  two  problems  associated  with  the  numerical  use  of  integral 
equations  in  scattering  calculations.  First,  the  kernels  of  Lippmann-Schwinger  equations 
are  dense,  except  when  the  background  is  extremely  attenuating;  since  iterative  techniques 
require  application  of  the  matrix  of  the  discretized  integral  operator  to  a  sequence  of  re¬ 
cursively  generated  vectors,  the  cost  of  the  procedure  is  prohibitive,  except  for  extremely 
small-scale  problems.  This  difficulty  was  overcome  almost  40  years  ago  via  the  observa¬ 
tion  that  the  free-space  Green’s  function  for  the  Helmholtz  equation  is  translation  invariant; 
appropriately  chosen  discretizations  of  Lippmann-Schwinger  equations  result  in  Toeplitz  ma¬ 
trices,  and  the  latter  can  be  rapidly  applied  to  arbitrary  vectors  via  the  FFT  (Fast  Fourier 
Transform),  resulting  in  algorithms  with  CPU  time  requirements  proportional  to  N  -log  (N), 
with  N  the  number  of  nodes  in  the  discretization  of  the  problem.  Various  forms  of  this 
approach  have  been  widely  used  in  electrical  engineering  and  other  holds,  under  the  name 
“k-space”  methods;  some  of  the  existing  codes  are  quite  fast,  even  for  discretizations  involv¬ 
ing  hundreds  of  millions  of  nodes.  However,  the  resulting  solvers  for  the  underlying  PDEs 
are  usually  not  very  accurate,  due  to  the  problem  discussed  in  the  following  paragraph. 

The  second  difficulty  associated  with  numerical  use  of  Lippmann-Schwinger  equations 
is  due  to  the  singular  character  of  the  Green’s  function  for  the  Helmholtz  equation;  in 
two  dimensions,  the  principal  term  of  the  singularity  is  of  the  form  log(r),  and  in  three 
dimensions,  it  is  of  the  form  1  jr.  As  a  result,  kernels  of  Lippmann-Schwinger  equations  are 
singular;  the  singularities  are  located  on  the  diagonal,  and  in  two  dimensions  are  of  the  form 

K(x,  y )  =  log (| or  -y |)  +  P(x,  y )  •  logflx  -  y |)  +  Q (x,  y),  (13) 

with  P,  Q  two  smooth  functions,  and  P(x,x)  =  0  for  all  x  G  R 2;  the  corresponding  form  in 
three  dimensions  is 

K(x,y)  =  t  1  ,+P(x,y)-  1  y  +  Q{x,  y).  (14) 

\x-y\  \x-y\ 

It  is  important  to  observe  that  in  most  cases,  we  do  not  have  access  to  each  of  the  functions 
P,Q  separately,  but  can  only  evaluate  the  whole  kernel  K  given  a  pair  of  points  (x,y).  In 
other  words,  standard  integration  techniques  (such  as  product  integration,  etc.)  can  not  be 
used  efficiently.  The  standard  procedure  in  the  literature  (referred  to  as  the  “singularity 
extraction”)  is  to  subtract  the  principal  singularity  and  treat  it  analytically,  and  apply  the 
trapezoidal  quadrature  rule  to  the  remaining  function.  Since  the  latter  is  not  smooth  (having 
infinite  derivatives  at  x  —  y),  the  procedure  converges  slowly,  normally  behaving  like  a  second 
order  scheme. 

We  introduce  a  class  of  quadrature  formulae  for  functions  of  the  form  (13)  in  two  di¬ 
mensions  and  (14)  in  three  dimensions.  The  approach  is  somewhat  related  to  the  Ewald 
summation  [6],  and  leads  to  quadratures  that  can  be  viewed  as  a  version  of  the  corrected 
trapezoidal  rule;  it  is  easily  combined  with  the  FFT  to  obtain  fast  algorithms.  While  in 
principle  corrections  of  arbitrarily  high  order  could  be  constructed,  in  practice  both  the 
complexity  of  derivation  and  the  number  of  corrections  grow  rapidly  with  the  order.  We 
have  designed  corrections  of  orders  4,  6,  8,  and  10;  they  require  1,  5,  13,  and  25  corrected 
nodes  respectively. 

The  paper  is  organized  as  follows.  In  Section  2,  we  summarize  several  well-known  math¬ 
ematical  facts  to  be  used  in  the  paper.  In  Section  3,  we  introduce  analytical  tools  to  be  used 
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in  the  construction  of  the  algorithm.  Section  4  describes  the  algorithm  in  detail,  and  contains 
a  complexity  analysis.  In  Section  5,  several  numerical  examples  are  used  to  illustrate  the 
performance  of  the  algorithm.  Finally,  Section  6  contains  generalizations  and  conclusions. 


2  Analytical  Preliminaries 

In  this  section,  we  summarize  several  well-known  mathematical  facts  to  be  used  in  the 
sections  below.  All  of  these  are  either  well  known  or  easily  derived  from  well-known  results. 


2.1  Notation 

For  an  integer  N  >  1,  the  two-dimensional  discrete  Fourier  transform  TN  is  a  mapping 
converting  a  two-dimensional  complex  sequence  a  =  {aJU2},  j\,j2  =  —TV,  into  another 

two-dimensional  complex  sequence  A  =  Ay,  k2  =  —TV, TV,  defined  by  the  formula 

N  N 

Aklk2=  Y  Y  ajlj2e-mk^e~mk^.  (15) 

ji=—N  j2=—N 

It  is  easily  verified  that  the  inverse  (JF^)-1  of  the  mapping  jFiV  is  given  by  the  formula 


( -F  )  {A)jij2  ~  ahh  ~ 


N 


N 


(27V  +  l)2 


Z^  Z^ 

ki=—N  k2=—N 


Aklk2e^^>kljl  e^^k2j2,  (16) 


with  j i  =  -N, ...,  N,  j2  =  -N, 

For  a  Helmholtz  equation 

W2(t)  +  k2(t)  =  0  (IT) 

in  two  dimensions,  the  potential  0  at  a  point  x  produced  by  a  unit  point  source  at  x0  is 
given  by  the  formula 

o(.r)  =  — |//(l(/.i|.r  -  .r0||),  (18) 

where  k  is  a  complex  number  such  that  Im(fc)  >  0,  and  Hq  is  the  Hankel  function  of  the  first 
kind  of  order  zero.  The  well-known  Sommerfcld  formula  states  that 

H0(kr)  =  -■  7  •  eiVW^x  •  eiXy  dX 

TT  oo  y/k2  -  A2 

for  any  k  G  C+,  r,  x,  y  >  0,  and  r  =  yj x2  +  y2  (see,  for  example,  [9]). 

Finally,  we  will  need  the  identity 

n  n  n— 1  j 

fj9i  =  fn  9k  -  5Z(/j+i  -  fj )  •  CY 9k )> 

j= 0  k= 0  j= 0  k= 0 


(19) 


(20) 


valid  for  two  arbitrary  finite  sequences  {fj},  j  =  0, 1,  2, ... ,  n„  {gj},  j  =  0, 1,  2, ... ,  n,.  By 
analogy  with  integration  by  parts,  (20)  is  normally  referred  to  as  summation  by  parts;  it  is 
easily  verified  by  a  substitution. 
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2.2  High-Order  Corrected  Trapezoidal  Quadrature  Rules  for  Sin¬ 
gular  Functions  in  One  Dimension 


For  a  function  /  :  [a,  b]  — >  R 1  and  integer  n  >  2,  the  n-point  trapezoidal  rule  Tn  is  defined 
by  the  formula 


(21) 


where 


h  = 


b  —  a 
n—  1  ’ 


(22) 


and  is  widely  used  as  an  approximation  to  the  integral  J^f(x)dx.  It  is  second  order  con¬ 
vergent,  as  long  as  the  second  derivative  of  /  is  continuous  on  [a,  b\.  In  other  words,  if 
/  G  C2[a,  b],  then 

f(x)dx  =  Tn{f)  +  0(h2).  (23) 


For  any  function  /  G  C2m+2[a,b ]  with  integer  m  >  1,  the  well-known  Euler- Maclaurin 
formula  (see,  for  example,  [1])  states  that  there  exists  a  real  number  £  with  a  <  C,  <  b,  such 
that 


/b  m  l  21  R  u2m+2  R 

f(x)dx  =  T„(/)  +  |jl2yf(/(2,-1,(6)  -  /<*-*>(«))  -  (2m+2"f2/|2m+2)(e), 


(24) 


where  Bk,k  =  0,1,2,...  denote  the  Bernoulli  numbers  (see  [1]).  It  is  easily  seen  from  (24) 
that  for  any  function  /  G  Cm[a  —  mh,  b  +  mh\  with  integer  m  >  3,  it  is  possible  to  construct 
quadratures  of  the  form 


m—  1 
2 

Tp™(f)  —  T„(f)  +  h  V  (/(&+  kh)-  f(a  +kh))f%,  (25) 


where  f3™  are  coefficients  such  that 

[  f(x)dx  =  !>"„(/)  +  0(/!-”+1),  (26) 

J  a 

with  h  defined  by  (22).  The  quadrature  T^m  is  referred  to  as  the  (m  +  l)</l-order  endpoint- 
corrected  trapezoidal  rule;  for  any  given  k  and  m,  where  m  >  3,  and  is  odd,  the  coefficient 
can  be  obtained  via  a  direct  calculation  (see  [7]). 

While  (21)  and  (25)  are  widely  used  for  the  numerical  integration  of  smooth  functions, 
their  use  for  singular  integrands  tends  to  encounter  difficulties  (see,  for  example,  [11]).  In 
[7],  a  class  of  quadrature  formulae  is  constructed  for  approximating  integrals  of  singular 
functions  of  the  form 

f(x)=(/>(x)s(x),  (27) 
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where  <p  G  Cm[—b  —  mh ,  b  +  mh]  with  m  >  3,  and  s  is  an  integrable  function  on  [—b,  b]  with  a 
singularity  at  zero.  For  integers  n  >  1,  p  >  1,  and  odd  m  >  3,  the  quadrature  UThpm  for  the 
functions  with  separable  singularities  (i.e.,  functions  defined  by  (27))  is  given  by  the  formula 


CW»(/)  =  r'JU/)  +  E  fVO'i)- 


(28) 


J=~P 


In  (28), 


=  h 


( 


E  fw 


J=-n 


f(~b)+f(b) 


\  j/o 


/ 


m—  1 
2 


+  ft  E  (f(<>  +  kh)-f(-b  +  kh))/3Z, 


k=-'~ 


(29) 

where  the  coefficients  f3™  can  be  obtained  via  direct  calculation  (see  [7]),  h  —  and  the 
coefficients  rj'  can  be  obtained  by  solving  the  system  of  linear  equations 


p  ,-b 

E  xTlri  =  /  (^*_1s(x))  dx  —  T'pm(x'l~1s), 

j=-P  J~b 


(30) 


with  Xj  =  jh ,  and  i  —  1,  2, . . . ,  2p  +  1.  The  quadrature  formula  (28)  is  of  order  2 p  +  2. 

Remark  2.1  While  (21)  is  the  standard  trapezoidal  rule,  and  (25)  is  the  trapezoidal  rule 
with  endpoints  corrections,  (28)  is  the  trapezoidal  rule  with  corrections  at  both  endpoints 
and  a  singular  point  inside  the  interval.  Here,  the  singular  point  is  the  center  point,  so  the 
scheme  is  sometimes  called  a  “center-corrected  trapezoidal  rule.”  The  coefficients  rj1  are 
called  correction  coefficients. 


Remark  2.2  The  only  difference  between  and  Tpm(f)  is  that  does  not 

contain  the  term  h  ■  /( 0),  which  may  become  infinite  for  /  of  the  form  (27). 


2.3  Toeplitz  Convolution 


This  section  introduces  two-dimensional  Toeplitz  convolutions  and  a  procedure  for  the  cal¬ 
culation  of  two-dimensional  Toeplitz  convolutions  via  the  two-dimensional  discrete  Fourier 
transform.  The  Toeplitz  convolution  a  *  b  of  finite  two-dimensional  complex  sequences 
a  =  {ajij2},  Ji, h  =  and  b  =  {bjlj2},  j1,j2  =  -2N,...,2N,  is  defined  by  the 

formula 


N  N 


(a  *  b)klk2  =  E  E 


a,- 


ljlj2bk1-j1,k2~j21 


ji=-N  j2=-N 


(31) 


where  k  \ ,  k2  =  The  well-known  convolution  theorem  states  that  the  Toeplitz 

convolution  a  *  b  is  equal  to  the  inverse  Fourier  transform  of  the  product  of  the  Fourier 
transform  of  a'  and  b ,  where  a'  is  a  two-dimensional  sequence  obtained  by  padding  the 
two-dimensional  sequence  a  with  zeros.  In  other  words, 


(o*%lfa  =  (rm)-'(r2N{a')-rm(b))kihi, 


(32) 


where  /ci,fc2  =  and  the  coefficients  of  the  two-dimensional  complex  sequence 

a'  =  {aii?:2}’ ?'i>  ^2  =  —  21V,  ...,21V  are  defined  by  the  formulae 

ai\i2  =  ahni  (33) 

when  —  N  <  ii,i2  <  IV,  otherwise 

=  0.  (34) 

Remark  2.3  While  direct  calculation  of  Toeplitz  convolution  (31)  leads  to  time  cost  of  order 
0(1V4),  which  is  prohibitive  for  large  scale  problems,  application  of  FFT  to  the  formula  (32) 
reduces  the  time  cost  to  0(N 2  ■  log  TV)  (see,  for  example,  [2]).  In  this  paper,  FFT  is  used  for 
the  fast  calculation  of  Toeplitz  convolution. 


3  Mathematical  Apparatus 

In  this  section,  we  introduce  analytical  tools  to  be  used  in  the  construction  of  the  algorithms. 

3.1  Endpoint  Corrected  Trapezoidal  Quadrature  Rules  in  Two  Di¬ 
mensions 

This  section  can  be  viewed  as  the  extension  of  results  of  Section  2.2  to  two  dimensions. 

For  a  function  /  :  [a,  b]  x  [a,  b]  — >  R 1  and  integer  n  >  2,  the  two-dimensional  trapezoidal 
rule  T2D  is  defined  by  the  formula 


n—  1  n—  1 

«)  =  EE  /(a  +  ih,  a  +  jh )  •  h2  ■  (35) 

i= 0  j= 0 


where 


h  = 


n  —  1  ’ 


(36) 


and  faj  equals  1  in  the  interior  of  the  square  [a,  b]  x  [a,  b] ,  equals  \  in  the  interior  of  the  edge, 
and  equals  t  on  the  corners  of  the  square. 

Further,  if  /  G  C2[a,b]  x  [a,b],  then 


[  [  f(x,y)dxdy  =  T™(f)+0(h2).  (37) 

J  a  J  a 

The  proof  consists  of  a  straightforward  application  of  one-dimensional  trapezoidal  rule  (21) 
to  both  directions  in  two  dimensions,  and  thus  omitted. 

The  following  Lemma  provides  the  two-dimensional  version  of  the  Euler-Maclaurin  for¬ 
mula. 
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Lemma  3.1  Suppose  that  a  function  f  G  C2m+2[a,  b]  x  [a,  b\,  integers  m  >  1,  n  >  2  .  Then, 


h2l+1  •  B2l  •  A 
(2/)! 


/.fr  i=n—  1  m 

/  /  /(x,j/)dxdj/  =  T^(/)  +  El 

l'a  i=0  «=1 

d2l~^  d2l~^  \ 

■^rif(b1a  +  ih)-^rif(a1a  +  ih)-^rTf(a  +  ih1b)-^rif(a  +  ih1a)j 


h2l+21'  ■  Bo,  ■  Bov 


1=1  l'= 1 


Q2l+2l'-2 

dx 2l~1  d  y2 


-£  E  -  (2/!)g,n^  ( (m *)+/(«> «)-/(«. «)) )  +o(/^2m+2), 


A  = 


where  T2D(f)  is  defined  in  (35),  h  is  defined  in  (36), 

1  0  <  i  <  n  —  1 

1/2  i  —  0  or  i  —  n  —  1  ’ 

and  Bk,  k  =  0,1,  2, ...  denote  the  Bernoulli  numbers. 

Proof.  The  Euler-Maclaurin  formula  (24)  can  be  rewritten  as 

h2lBo 


(38) 

(39) 


f(x)  dx  =  TM)  +  E  ~  f(2‘-VIM)  +  0(ft2m+2).  (40) 


Hence, 

f(x,y)  dxdy  = 

rb  n_1  m  U21  to 

/  (£/(«  +  !fc,9) ■  ft  +  E  -h7#(/(2'_1,(f>)  -  /(2|-1)  M)  +  0(h2m+2))  dy.  (41) 

The  conclusion  of  the  Lemma  follows  immediately  from  applying  the  formula  (40)  to  the 
integrals  in  (41).  □ 


The  following  Lemma  provides  a  (2 m  +  2)^l-order  endpoint-corrected  trapezoidal  rule  in 
two  dimensions. 

Lemma  3.2  Suppose  that  a  function  f  G  C2m+2[a,b]  x  [a,b]  with  integers  m  >  1,  n  >  2, 
h  =  ^  .  Then, 

n—  1  7 

f  f  f(x,y)  dxdy  =  T^{f)  +  0(h2m+2),  (42) 

J  a  J  a 

In  (f2), 


rri2,T)  ftX 
Q2m+1 


(/) 


n—  1  m 

r,f  (/)  +  V  E  E 


i=0  k=—m 


(ft2m+i 


A- 


m  m 

(/(6  +  Jfe/i,*/i)  +  /(*/i,6  +  A:/i)^/(a  +  A:/i,*/i)-/(*/i,a+A:/i)))  +h2  ^  ^  (^m+1^,m+1- 

k=—m  k'=—m 

(/(a  +  kh ,  a  +  A//r)  +  /(&  +  kh ,  6  +  k'h )  —  /(a  +  kh,  b  +  A//i)  —  /(6  +  kh,  a  +  A/h))  j ,  (43) 
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T2D ,  (3lm+1,  and  (3i  are  defined  by  (35),  (25)  and  (39). 

Proof.  It  is  easily  observed  from  the  formulae  (24),  (25)  that 

h2l+1B ■ 


S  en> 


SL  f(2l~ i)(n\  _ 


a  = 


]T  f(a+kh)./d2km+1-h2. 

k=—m 


Therefore, 


and 


h 2  ■  f(a  +  kh,  a  +  ih)  ■  (3l 


2m +1 
k 


k=—m 


^  h*+'B2l  „  ,, 


EE 

/=1  l'= 1 

m  l2J'  r 
'l  -°2 11 


h2l+2V  ■  b21  ■  b21,  (  d2l+2l'~2  f(  ^ 
'!)  vaxa-15j/a'-l/^a7 


(2/!)(2Z 

m 

E 


k=—m 


dy 


m  ,  m  U2l'~  1  R  02i'-l 

E^r1- 


(2/d)  fy2* 

/i2  •  /3lm+1/3l,m+1/(a  +  kh,  a  +  k'h) 


k=—m  l 

m  m 


k=—m  k'=—m 

Now,  (43)  follows  immediately  from  the  combination  of  (45),  (46). 


(44) 


(45) 


(46) 

□ 


3.2  High-Order  Center  Corrected  Trapezoidal  Quadrature  Rules 
for  the  Singular  Functions  in  Two  Dimensions 

The  following  Lemma  provides  an  estimate  of  the  difference  between  the  integral  and  the 
end-point  corrected  trapezoidal  quadrature  for  functions  of  the  form  x2p+2  ■  log(x2  +  y2). 

Lemma  3.3  Suppose  that  n  is  a  positive  integer,  a,  h  are  two  positive  real  numbers  such 
that  h  =  a/n,  and  integers  m,  p  are  such  that  m  >  p  +  1  >  1.  Then, 

[  f  (: x2p+ 2  ■  log(a;2  +  y2))  dx  dy  =  T^S+i  {x2p+2  ■  log(a;2  +  y2))  +  0(h2p+i ),  (47) 

J  —a  J  —a 

with  i  defined  by  (43). 

Proof.  Here,  we  prove  the  case  of  p  =  0  as  an  illustration;  the  proof  for  p  >  0  is  quite 
similar.  For  simplicity,  we  will  be  denoting  x2  ■  log(x2  +  y 2)  by  g(x,y).  Then, 

dg(x,y)  =  2 x2y 

dy  x2  +  y2’  1  ’ 
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and 


d3g(x,y )  16 x2y3 


dy3 


12  x2y 
( x2  +  y 2)3  (a;2+y2) 


2’ 


d6g(x,y)  7680 x2ye  |  1 1520a;  V  4320a;V  |  240a;2 

”1”  t  ^  t  ^  rrr  + 


ch/6 


(X2  +  !/2)6  (l2  +  !/2)5  (l2  +  !/2)4  (X2  +  V2)3' 


(49) 


(50) 


Replacing  the  derivatives  in  (24)  with  the  appropriate  finite  differences,  we  rewrite  (24)  in 
the  form 


J  g(x )  dx  =  Tn(g )  +  —  h  •  (■ -g(b  +  h)  +  ^(a  +  h)  +  g(b  -  h)  -  g(a  -  h))  + 

Cl  •  fc4  •  (j<3>(6)  -  g<3>(a))  -  C2  ■  tis  ■  g<6>«) 

=  T$,(g)  +  O  •  h‘  ■  (g<3>(6)  -  9,3|W)  +  C2  •  A6  •  J(6)K),  (51) 

where  C\  and  C2  are  two  constants  independent  of  g,  and  £  G  [a,b\.  Therefore,  for  any 
fixed  x  and  g(x,y)  =  x2  ■  log(a;2  +  y2),  the  error  e(x)  of  the  end-point  corrected  trapezoidal 
quadrature  in  the  y  direction  is 


16a;  a 


2^3 


12a;2a 


<x)  =  j_a9(x,y)  dy  -  Tft(g)  =  C \  •  h4(^2  +  a2)3  -  (a;2  +  a2)2j  +  C2  •  h 6  •  g{6\x,0-  (52) 
Summing  up  all  the  errors  along  the  x  axis,  we  obtain 

pa  pa  n 

(. x 2  ■  log(a;2  +  y2))dxdy  —  T23D,n(x2  ■  log(a;2  +  y2))  ~  2  ■  ^  h  ■  e(ih ) 
16(ih)2a3 


e  = 


'  —a  J  —a 
n 


i= 1 


<  h  •  ^  (  C4  •  h 


1=1 


((ih)2  +  a2)3  ((* 

,  „  ,4  /'a  /  16a;2a3 

^  Ol  ll 


^_)+ft.gC2.A».  1^(01, 


12a;2  a 


o  y  (a;2  T  a2)3  (a;2  +  a2)2 


da;  +  C2  •  h7  — 

Z— /  (  n 


240 


C*!  •  (3  —  7r)  •  /a4  T  C*2  ■  /r6  ■  log(h)  ~  0(h4). 


(53) 

□ 


Theorem  3.4  below  is  an  extension  of  the  formula  (28)  to  two  dimensions.  For  an  integer 
p  >  0,  it  supplies  a  (2 p  +  4)<h-order  center-corrected  quadrature  formula  on  R2  for  the 
functions  of  the  form 

f(x,y)  =  <t>(x,y)-s(x,y),  (54) 

where  (j)  :  K2  — >  R,  and 

s(x,  y)  =  log(a;2  +  y2)  +  P(x,  y)  •  log(a;2  +  y2)  +  Q(x,  y),  (55) 

with  P,  Q  two  smooth  functions,  and  P(0,  0)  =  0.  Suppose  that  n,  m  are  positive  integers, 
and  a,  h  are  two  positive  real  numbers  such  that  h  =  a/n.  We  define  T'2^^+1{f)  by  the 
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formula 


r'J2»”+1(/)=  5]  f(ih,jh)-h2-/3i:i  +  h2Yl  'll 


i=—n  k=—m 


m  m 

(/(a  +  kh,  ih)  +  f(ih,  a  +  fc/i)  —  /(— a  +  kh,  ih)  —  f(ih,  —  a  +  kh))  j  +h2  Y^  ^/3^m+1/3^Jr 

k=—m  k'=—m 

(/(—a  +  kh,  —a  +  A//i)  +  /(a  +  kh,  a  +  k'h )  —  /(—a  +  kh,  a  +  k'h)  —  f(a  +  kh,  —a  +  k'h ))  j , 


where 

M  =  {i,j  e  Z  :  |i|  <  n, \j\  <  n,  (■ i,j )  ^  (0,0)},  (57) 

A- {1/2  M-»,1  (58) 

in  (56),  coefficients  equal  1  in  the  interior  of  the  square  [—a,  a]  x  [—a,  a],  equal  |  in  the 
interior  of  the  edge,  and  equal  }  on  the  corners  of  the  square,  and  /?2m+1  are  defined  in  (25). 

Theorem  3.4  Suppose  that  d>  e  C2p+2(R2)  withp  >  0,  integer  m  >  p+ 1,  and  s  is  a  singular 
function  on  R2  with  a  logarithmic  singularity  at  (0,0),  he.,  of  the  form  (55).  Suppose  further 
that 

U2^2m+fifi-s)=T,2°mn+1(cj)-s)  +  Y  r^iihjh),  (59) 

where  ■  s )  is  defined  by  the  formula  (56), 

W  =  (i,  j  £  Z  :  |i  +  j|  <  p  and  \i  —  j\  <  p},  (60) 

and  the  coefficients  t(-  in  (59)  satisfy  the  system  of  linear  equations 

_  [  [  (xl'~1y3'-1s{x,y))  dxdy  -  y3' s) ,  (61) 


Y  xi  ly\ 

( i,j)ew 


' .  T-  ■  = 
3  l3 


with  Xi  =  ih,  ijj  =  jh,  and  ( i',j' )  e  H ,  where  H  =  {i' ,j' ,  &  Z  :  i'  >  1, /  >  1,*'+/  <  2p  +  2}. 
Then, 

f  f  (fi(x,y)  ■  s(x,y))  dxdy  =  •  s)  +  0(h2p+4).  (62) 

J  —a  J  —a 

Proof.  Applying  the  Taylor  expansion  to  the  function  fi(x,  y)  at  the  point  (0,  0)  we  have 


fi(x,y)  =  P(x,y)  +  R(x,y), 


where 


2p+l  j 


P(x,v)  = 

3=0  i=0  V 


xiyj  idxidyj-*<l>(X'y^x=0’v=0' 
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and 


R(x,y)  = 


— - £  (^P  +  2  ]xi„2P+2-, 


(2p  +  2)!  y 


i= 0 


^2p+2 


(65) 


'  —a  d  —a 


< 


’  —a  J  —a 


'  —a  J  —a 


where  £i,G  £  [—a,  a]  x  [—a,  a].  Thus, 

pa  pa 

1/  /  (4>(x,y)  ■  s(x,y))dxdy -U^m+iicp- s)\ 

(P(x,  y)  ■  s(x,  y))  dx  dy  -  f/^2m+i  (P  ■  s)  |  + 

(. R(x,y )  ■  s(x,y ))  dxdy  -  C/^2m+i(-R  •  s)|. 

I 

Now,  we  make  the  following  three  observations.  Due  to  (61), 

[  [  (P(x,y)  ■  s(x,y))dxdy-Ulhpm+1{P  ■  s)  =  0, 

J  —a  J  —a 

and  due  to  (59), 

pa  pa 

\  /  (R(x,y)- s(x,y))dxdy-U?h02m+i(R-s)\ 


(66) 


(67) 


1  —a  J  —a 


<\  [  [  (R{x,y)-s(x,y))dxdy-T'p£.i(R-s)  |  +  |  ^  yg  R(ih,  jh)\. 
J~a  J ~a  ( i,j)ew 

Due  to  formulae  (47)  and  (61),  it  is  obvious  that 

I  /  f  (R(x,y)  ■  s(x,y))dxdy -T,2^+i(R- s)\  ~  0(h2p+li), 

J  —a  J —a 

and  the  coefficients  r/J  are  of  the  order  h.2,  hence, 


(68) 


^  T!‘1R(ih,jh)\~0(k2’'+i 
(i,j)ew 


(69) 


(70) 


Finally,  the  conclusion  of  the  Lemma  follows  from  the  combination  of  (66)  -  (70). 


□ 


3.3  High-Order  Center  Corrected  Trapezoidal  Quadrature  Rules 
for  the  Green’s  Function  for  the  Helmholtz  equation 

In  this  section,  we  prove  Theorem  3.10,  which  is  the  principal  analytical  tool  of  this  paper. 
Theorem  3.10  describes  the  10^' -order  center-corrected  trapezoidal  quadrature  formulae  for 
the  Hankel  function.  It  can  be  viewed  as  a  special  case  of  Theorem  3.4  with  p  =  3,  and  the 
singular  function  s  the  Hankel  function  Hq. 
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In  the  remainder  of  this  paper,  we  will  be  using  the  following  notation.  For  any  k  E  C+ 
and  h  >  0,  we  will  define  the  complex  numbers  D0,  D±,  D2 ,  D3,  D 4,  D5,  via  the  formulae 


f*  oo  poo 


Dn  = 


'  — oo  J  — OO 


H^\kr)  dx  dy  —  ^  H^\k\J(ph)2  +  (qh)2)  ■  h 2,  (71) 

(p, 9)^(0, o) 


f*  oo  /*oo 


Do  = 


'  — oo  — oo 


r»oo  /»oo 


'  — oo  J  — oo 


H^\kr)x2  dx  dy  —  ^  H^\ky/ ( ph )2  +  (gh)2)  •  (ph)2  •  h2,  (72) 

(P, 9)^(0, 0) 

H^\kr)x4  dx  dy  —  ^  H^\koJ (ph)2  +  (gh)2)  •  (ph)4  ■  h2,  (73) 

(p,  9)^(0,' o) 


r»00  POO 


Do  = 


'  — OO  d  — oo 


H^\kr)x2y2  dx  dy  —  ^  H^\ky/ (ph)2  +  (qh)2)  ■  ( ph)2(qh )2  •  h2,  (74) 

(p, 9)^(0, o) 


poo  /»oo 


F>4  = 


r  — OO  d  — oo 


H^\kr)x6  dx  dy  —  ^  H^\k^/ (ph)2  +  (qh)2)  ■  (ph)6  ■  h 2,  (75) 

0, 9)^(0, 0) 


r»00  /»00 


As  = 


H^\kr)x4y2  dx  dy  —  ^  H^\ky/ (ph)2  +  (gh)2)  •  ( ph)4(qh )2  ■  h2 .  (76) 

(P,  9)^(0, 0) 

The  following  Lemma  is  a  simple  consequence  of  the  Sommcrfeld  formula  (19). 


'  —  OO  J  — oo 


Lemma  3.5  For  any  k  E  C  ,  r,  x,  y  >  0,  and  r  =  J  x2  +  y 


1  r°° 

H0(kr)  =  - 


.  gi-^-(Vfc2-A2 -A)-x  .  ei-^-(Vfc2-A2+A)-y 


(77) 


J-oo  V*:2  -  A2 
where  r 2  =  x2  +  y2,  x,y  >  0. 

The  following  two  technical  lemmas  follow  immediately  from  the  Sommcrfeld  formula  (19). 
Lemma  3.6  For  any  k  E  C+,  and  a  >  0, 

4  r  1 


H0(kr )  dxdy  = 


7T  J- oo  aA2  -  A2 

ei-^-(Vfe2-A2-A)-a  _  ■]_  g«-^-(yfc2-A2+A)-a  _ 


r  —a  «/  —a 
V2 


i  ■  &  ■  (Vk2  -  A2  -  A)  i  •  ^  •  (Vh2  -  A2  +  A) 


dA, 


(78) 


with  r  =  \J x2  +  y2. 

Proof.  Substituting  (77)  into  the  left  side  of  (78),  and  changing  the  order  of  integration, 
we  obtain 


'  —a  d  —a 


H0(kr )  dxdy  =  4 
4  f00  dA 


o  jo 


7T  d-oo  Vh2  -  A2 

4  r00  i 


H0(kr )  da;  dp 
f  ei-&-Wk*~ Aa-A)-*  dx. 


g*-A^-(\/fc2— A2— A)-a  gi-^-(v/fc2-A2+A)-a  


ei-4-WW~  *+Dydy 


V2, 


71 


l-oo  Vk2  -  A2  i.f  (v^l-A)  i  •  ^  •  (Vh2  -  A2  +  A) 


dA.  (79) 
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□ 


Lemma  3.7  For  any  k  G  C+,  integer  n  >  1,  and  a  >  0 

n  n 

E  E  H0(ky/[ph )2  +  ( qh )2)  •  Ii2  ■  (3, 

1 


-'pq 


(80) 


p=—n q=—n 

4  f00 


•M- 


ri-F?L-(\/k2—X2—X)-a  i  i  -i  _ 

e  2  _  ~  1  _  i  ,  ± 

ei-^  (VkZ=>?-X)  h  _i  2  2 


7T  .Moo  \/&2  -  A2 

X-^S--(\/k2— A2+AVa  i  i  i  _ 

- _E  _  I  +  le‘-4  (Vi^+A).  rfA 

^i-2^  ■  (\/  k'2  —  X2 +X)-h  _  2  2  2  ' 


(81) 


with 


h  =  a/n ,  (82) 

and  Ppq  equals  1  m  t/je  interior  of  the  (2 n  +  1)  x  (2n  +  1)  square,  equals  \  in  the  interior  of 
the  edge,  mid  equals  |  on  the  corners  of  the  square. 

Proof.  The  trapezoidal  sum  (84)  over  the  domain  [—a,  a]  x  [—a,  a]  is  equal  to  four  times 
the  trapezoidal  sum  over  the  domain  [0,  a]  x  [0,  a].  In  other  words, 


n  n 


Ho\kV(ph)2  +  (qh)2)  ■  h2  ■  [3pq 


p=—n q=—n 

n  n 


4  '  E  E  H«\W(ph)2  +  («A)2)  •  h2  ■  I3'm, 

p= 0  q= 0 


(83) 


where  fi'pq  equals  1  in  the  interior  of  the  (n+ 1)  x  (n+ 1)  square,  equals  |  in  the  interior  of  the 
edge,  and  equals  |  on  the  corners  of  the  square.  Substituting  (77)  into  (83),  and  exchanging 
the  order  of  integration  and  summation,  we  obtain 


E  E  H0(k v7 iph )2  +  (qh)2)  •  h2  ■  (3pq 

p=—n q=—n 

4  dX 


(84) 


vr  \Jk2  -  A2 


•"ME 


—  1  I  „i- ^  ■  (V  k2  ~  X2  ~  A)- 

i ■  E  .(y/k2  —  X2—X)-p-h  1  e  _ 


e  2 


p=0 


n  1  ,  X-^-(\/k'2-X2+X)-a 

j.E.(v^2XA2+A).q./l  _  l  +  e  2 _ 


Se*'  2 
<?=0 


(85) 


vr  M/c2  -  A2 


•  /T 


r.i-2^-(\/k,‘2—X2—X)-a  i  i  -|  _ 

e  2  1  j  ~  1  _  I  ,  iei-^.(MG3A2-A).( 

ei.^.(Vk^x^-\)-h  _  i  2  2 


ei.E.(Vfc2— A2+A)-a  _  x  x  X  ^ 


_  _  +  _ei-^-Wk2-X2+X)-a  ]  dX 
x^-Wk2-X2+X)-h  _  i  2  2 


(86) 

□ 
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Remark  3.1  As  a  — >  oo,  the  exponential  terms  e^Y'^k'2  X2±F-a)  in  (78)  and  (81)  tend  to 
zero;  this  fact  will  be  used  in  Lemma  3.8  below. 

The  following  lemma  provides  an  analytical  form  for  the  difference  between  the  inte¬ 
gral  (78)  and  the  trapezoidal  sum  (81).  Its  proof  consists  of  combining  Remark  3.1  with 
(78),  (81). 

Lemma  3.8  For  any  k  G  C+  and  h  >  0, 


f*  oo  r  oo 


oo  oo 


H^\kr)dxdy-  Y  Y  n£\ky/(ph)2  +  (qh,)2)  ■  h2 


'  — oo  J  — oo 


/■°°  d\ 

(  1 

/l2  e*aift  +  leia2/!  + 

/-oo  Vk2  -  X2 

V(mi)(m2) 

4  giai/i  _  7  giazh  _  7  I 

where 


p=— oo  q—— oo 

4  i  ua  /  -l  it  c-~r±c-~r±\  (QY) 

(88) 

2  \ . 2  ' . /■  (89) 

The  following  Lemma  follows  immediately  from  Lemma  3.8.  It  supplies  an  analytical 
form  for  the  difference 


r  =  \J  x1  +  y2, 


a !  =  -  A),  a2  =  +  A). 


f*  oo  /»oo 


'  — oo  — oo 


1yJ  1  H${k\/x2  +  y2)  j  da;  dy—  Y^  (iP^T  1  H{)(k\f  ( ph )2  +  (g/r)2)  j  -/?2, 

(90) 


(P, 9)^(0, 0) 


with  i  =  l,j  =  l. 


Remark  3.2  (90)  is  the  right  hand  side  of  the  equation  (61),  and  thus  is  directly  used  in 
the  calculation  of  coefficients  rh.  Direct  numerical  subtraction  of  the  integral  and  the  sum  in 
(90)  leads  to  loss  of  accuracy  because  of  cancellation  errors,  especially  when  i ,  j  are  relatively 
large.  Lemma  3.9  below  and  Lemmas  6.1  -  6.5  in  Appendix  A,  provide  analytical  formulae 
for  (90)  with  (i,j)  =  {(1, 1),  (3, 1),  (5, 1),  (3,3),  (7, 1),  (5,3)},  i.e.,  D0  -  D5  defined  by  (71)  - 
(76),  so  that  cancellation  errors  are  reduced. 

Lemma  3.9  For  any  k  G  C+  mid  h  >  0, 

D0  = 


r»QO  /»QO 


'  — OO  J  — OO 


H^\kr)  dxdy-  Y  Ho\W(ph)2+  (qh)2)  ■  h2 

(p,<2)tR  0.0) 


d\ 


h2 


Di<y.\h  _j_  gia.2h 


vr  oo  Vk2  -  A2  V  (*ai)  (*a2)  2  (eia^h  -  1) (eia*h  -  1) , 

4  d,\  (l — |j)h2  +  zi  +  z2  +  i\ j\\h{y\  —  y2)  +  ViV2 


7T  oo  y/k2  -  A2  (iai)(ia2)(l  +  £i)(l  +  x2) 

where  the  complex  numbers  x i ,  x2,  yi,  y2,  24,  z2  are  defined  by  the  formulae 

eiaih  _  ^ 


Xi  = 


ioiih 


-i  =  E 


(iaq/i)71 
^  (n+  1)!’ 


eiQ2/l  -  1  ^  (ia2/r)r 

^2  =  — - - 1  =  2^ 


fa2/?. 


^  (n  +  1)!’ 


(91) 

(92) 

(93) 
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y  i  =  xi 


zi  =  yi- 


ia\h  ( ia\h)r 

2~~  ~  ^ 

(iaih) 


n= 2 

2  00 


^  (w  +  l)!’ 


6 


E 


£<b+1>1’ 


1/2  =  X2 


Z2  =  1/2  - 


ia2h 


2  E 


(ia2h)n 
nS  ("+1)1’ 


(ia2  hfi 


E 


(■ ia2hy 


,t=s(«+d!’ 


(94) 

(95) 


and  r,  on,  and  a2  are  defined  by  (88),  (89). 

Proof.  Substituting  (77)  into 

Y  H^\ky/{ph)2  +  ( qh )2)  •  h2, 

(p>«)^(0,0) 

and  exchanging  the  order  of  integration  and  summation,  we  obtain 

Y  ( ph )2  +  (qh)2)  ■  h 2 

(p,?)^(o,o) 

A  f°°  dX  fh2eiaih  +  leia2h  +  l  h2 
~  7r  ^/fc2  -  A2  V  4  -  1  e*“2/l  -  1  4 

_  A  d\  h2  eiaih  +  eia2h 

”  vr  J V&2  -  A2  '  T  '  (e“1'“  -  l)(eiQ2h  -  1) ' 


(96) 


(97) 


Now,  (91)  follows  immediately  from  the  combination  of  (97),  (78),  and  Remark  3.1.  Substi¬ 
tuting  (93)  into  (91),  we  obtain 


h 2 


Dicx.\h  _|_  ^.OL^h 


[iai){ia2) 

1 


2  (giaifc  _  l)(eia2/i  _  1) 

1  2  +  iot\h{\  T  x i)  T  ia2h(l  T  xfi) 
(iafi)(ia2)  2  (ia1)(ia2)(  1  +  xi)(l  +  x2) 

Xi  +  x2  +  Xia;2  —  ifcKih(l  +  xfi)  —  \ia.\h(l  +  xfi 
{iafifiafifi  +  Xi)(l  +  x2) 

Finally,  (92)  follows  from  the  combination  of  (93),  (94),  (95)  and  (98). 


(98) 

□ 


Remark  3.3  Introducing  the  notation  z  =  we  rewrite  D$  in  the  form 
H^\kr)  dxdy-  Y  ^(^vW  +  W)  ■  h2 

(p>5)t^(0,0) 

4/r2  f00  dz  f  1  1  ei^-{Vl-z2-z)kh  _J_  ei^-(\/l-z2+z)kh 

VT  ,/_oo  \/l  —  ^2  V  (^2  —  0.5)  •  (kh)2  2  ^ei^(Vl-z2-z)kh  _  ^^eififi(y/l-z2+z)kh  _  ^ 

(99) 

Thus,  .D0  is  entirely  determined  by  k  and  h,  and  is  of  the  form  h2  ■  f{k  ■  h).  Similarly,  Di  is 
of  the  form  h 4  ■  f(k  ■  h);  D2  and  D3  are  of  the  form  h 6  •  f(k  ■  h );  H4  and  D5  are  of  the  form 
hs  ■  f(k  ■  h ).  In  other  words,  except  for  the  multiplicative  factors  (h2,hr,  h 6,  or  /i8 ) ,  Dq  -  D$ 
only  depend  on  the  product  k  ■  h. 
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Remark  3.4  Even  when  Lemmas  6.1  -  6.5  are  used,  a  certain  loss  of  accuracy  in  the  cal¬ 
culation  of  Di  -  D, 5  is  encountered  (see  Remark  3.2  above).  Thus,  evaluating  D0  in  double 
precision,  one  obtains  roughly  13  digits;  for  Di  one  gets  9  digits,  and  D2  D3  DA  D5  yield 
even  fewer  digits. 

To  avoid  this  difficulty,  we  utilized  extended  (real  *32)  precision  to  precompute  the 
coefficients  D0  -  D5  for  values  of  kh  at  appropriately  chosen  nodes  on  the  boundary  of  the 
square  11  =  [0, 1]  x  [0, 1]  in  the  complex  plane,  and  used  Lagrange  interpolation  to  evaluate  D0 
-  D\ 5  for  arbitrary  points  in  11  to  13  digits  (see  [8]  for  a  detailed  description  of  the  technique). 
Thus,  in  all  of  our  numerical  experiments  reported  in  Section  5  below,  the  coefficients  D0  - 
D: 5  were  obtained  by  interpolation,  rather  than  computed  “from  scratch” . 

Now,  we  are  ready  to  formulate  Theorem  3.10,  which  is  the  principal  analytical  tool  of 
this  paper  (together  with  Lemmas  6.6,  6.7,  6.8).  Theorem  3.10  describes  the  10*ft-order 
center-corrected  quadrature  formula  for  the  Green’s  function  for  the  Helmholtz  equation  in 
two  dimensions;  this  theorem  is  a  special  case  of  the  high-order  center-corrected  trapezoidal 
rule  for  singular  functions  in  two  dimensions  (see  Theorem  3.4)  with  p  —  3,  and  s(x,y)  = 
//{■ l>  (k \J x'2  +  y2) .  The  4th-order,  6^-order,  and  8*h-order  center-corrected  quadratures  are 
similar  and  listed  in  Appendix  B  (see  Lemma  6.6,  6.7  and  6.8).  All  the  proofs  are  quite 
similar  to  that  of  Theorem  3.4,  and  are  omitted. 


Theorem  3.10  Suppose  that  n  >  1  is  an  integer,  and  a,  h  are  two  positive  real  numbers 
such  that  h  =  a/n.  Suppose  further  that  0  :  R2  — »  C  is  a  function  such  that  f  e  c8  (R  x  R) , 
and  that  f  is  zero  outside  the  square  [—a,  a]  x  [—a,  a].  Then,  for  any  k  G  C+, 


[  1  <j>(x,y)  •  H0(ky/x2  +  y2)  dxdy  =  •  H0)  +  0(h 10). 

J  —a  J  —a 

(100) 

In  (100), 

UA<t>  ■  Ho)  =  T'2D(f  -H0)+Y,  Tpq&iph,  qh ), 

(101) 

p,q€S 

where 

S  —  {Pi  Q  €  Z  :  \p  +  q \  <  3  and  \p  —  q\  <  3}, 

(102) 

T/2D(0  -H0)=  {^Ph,  Qh)  •  ffftfcvW  +  W))  •  h 2, 

(103) 

(P, 9)^(0, 0) 

and 

h  t~\  49  D\  7  D2  3  D3  1  Di  1 D3 

ffio  -  0  _  187A  +  9~hd  +  2"/H  _  187^  _  27A  ’ 

(104) 

3  D\  13  D2  19  D3  1  Di  7  D§ 

T±io  ~  ro±i  “  4  h2  48  h4  24  hA  +  48  h 6  +  24  /r6  * 

(105) 

3  D\  1  D2  1  D3  1  Di  1  D§ 

r±20  -  ro±2  -  40  h2  +  12  hi  +  24  hi  120  he  24  h6  , 

(106) 

_  1  Dl  1-^2  1  Di 

±3°  0±3  180  h 2  144  hA  720  h 6 ’ 

(107) 

_  5  D3  1 D5 

T±1±1  ~  12  hA  6  he  ’ 

(108) 
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Figure  1:  The  25  correction  nodes 


h  _  h  _ 

r±l±2  ~  T±2±l  ~ 


48  hA  +  48  h6  ' 


(109) 


Remark  3.5  For  simplicity,  we  assume  here  the  function  0  to  be  zero  outside  the  square 
[—a,  a]  x  [—a,  a].  Thus,  the  endpoint  corrected  trapezoidal  rule  Tl2^m+i  (4>-H0)  in  Theorem  3.4 
reduces  to  the  standard  trapezoidal  rule  T/2D(0  •  H0),  and  the  integral  and  the  sum  on  the 
square  [—a,  a]  x  [—a,  a]  are  identical  to  those  in  R 2.  This  simplification  allows  the  direct 
use  of  the  analytical  formulae  for  Do  -  D$  (see  Lemma  3.9  above  and  Lemmas  6.1  -  6.5  in 
Appendix  A). 


Remark  3.6  Combining  Remark  3.3  with  the  definitions  (104)-(109),  we  observe  that  each 
of  the  coefficients  r^q  in  (104)-(109)  has  the  form  h2  ■  f(k-  h);  we  will  refer  to  the  coefficients 
Tpq  as  correction  coefficients. 

Remark  3.7  The  set  S  defined  in  (102)  contains  25  pairs  of  integers  ( p,q );  in  other  words, 
corrections  at  25  points  around  the  singularity  are  required  to  construct  a  10^-order  quadra¬ 
ture  formula  (see  Figure  1).  In  general,  for  any  integer  p  >  0,  2 p2  +  2p  +  1  correction  nodes 
are  needed  to  obtain  a  quadrature  of  order  2 p  +  4. 


3.4  Fast  Numerical  Application  of  Discretized  Lippmann-Schwinger 
Operators 

In  this  section,  we  combine  the  10t/l-order  quadrature  formula  for  the  integral  (100)  with  the 
FFT  to  obtain  a  fast  procedure  for  the  application  of  discretizations  of  the  operator  (12). 
We  will  denote  by  D  the  square  [—a,  a]  x  [—a,  a]  in  R2. 

Suppose  that  N  >  1  is  an  integer,  h  is  a  positive  real  number,  and  S'  is  a  set  defined 
in  (102).  Suppose  further  that  the  coefficients  are  defined  in  (104)-(109).  Then,  we 
define  a  two-dimensional  complex  sequence  H  =  {Hili2},  A,i2  =  —2 At, ...,  2N,  as  follows: 

Hni2  =  HoiW^h)2  +  M)2)  +  r^2/h2,  (110) 
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when  (*i,i2)  G  S'  and  (ii,i2)  7^  (0,0); 

i/oo  =  4/h2,  (111) 

and 

Hni2  =  tfo(*VM)2  +  M)2),  (112) 

otherwise.  We  define  a  complex  sequence  <h  =  {<hj1i2},  *i,*2  =  to  be  the  two- 

dimensional  sequence  defined  by  the  formula 

$*!«  =  (113) 

where  0  :  R2  — >  C  is  a  two-dimensional  c8-function  which  is  zero  outside  D. 

Lemma  3.11  Suppose  that  the  integers  n,  l±,  I2  are  such  that  n  >  l,  —n  <  h  <  n,  —n  <  l2  < 
n,  and  that  the  real  numbers  a,x,y  are  such  that  a  >  0,  —a  <  x  <  a,  —a  <  y  <  a.  Suppose 
further  that  0  :  R2  — >  C  is  a  c8 -function  which  is  zero  outside  the  square  [—a,  a]  x  [—a,  a]. 
Then  for  any  k  G  C+, 

[  f  4>{x',  y')  ■  H$\kyf{x  -  x')2  +  (y  -  y')2)  dx  dy 

J  —a  J  —a 

=  E  E  m-H(h-„m-n)  +  0(hla),  (114) 

—n<ii<n  —n<i2<n 


where  h  =  a/n,  x  =  l\h,  y  =  l2h,  the  two-dimensional  sequence  ii,i2  =  —n,  ...,  n 

is  defined  in  (113),  and  the  two-dimensional  sequence  H  =  { Hj u-2},  j\,j2  =  —2n,...,2n  is 
defined  in  (110)  -  (112). 

Proof.  Due  to  (100),  (101), 

[  <f>(x',  y')  ■  {kyj (x  -  x')2  +  (y-  y')2)  dx'  dy' 

J  —a 

=  ^  (, b[i\h,i2h )  •  H£\ky/ [fill  -  iih)2  +  (l2h  -  i2h)2)  ■  h2 
+  ^2  Tpq(f>(lih  +  pli,l2h  +  qh)  +0(h10),  (115) 

(p,q)£S 

where 

I’  =  {*i, i2  G  Z  :  |*i|  <  n,  \i2\  <  n,  (ii,i2)  ^  (116) 

and  S  is  defined  in  (102).  Now,  (114)  follows  immediately  from  the  combination  of  (115) 
and  the  definitions  in  (110)  -  (113).  □ 


Remark  3.8  Obviously,  (114)  is  the  Toeplitz  convolution  of  the  two-dimensional  sequences 
<h,  H,  and  as  such,  it  can  be  rapidly  calculated  via  the  FFT  (see  Section  2.3  above).  Thus, 

—n<i\<n  —n<i2<n 
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where  —n  <  l\  <  n,  —n  <  l2  <  n,  and  the  two-dimensional  sequence  =  {$h},  i,j  = 
—2 n, ...,  2 n,  is  dehned  by 


/  if  1*1  <  n  and  \j\  <  n 
\  0  if  |*|  >  n  or  \j\  >  n 


(118) 


Remark  3.9  For  any  point  x  outside  the  square  [—a,  a]  x  [—a,  a],  integral  (12)  is  approx¬ 
imated  via  the  standard  trapezoidal  rule.  This  approximation  is  10t/l-order  convergent,  as 
long  as  0  G  cw(R2)  . 


4  Description  of  the  Procedure 

This  section  describes  the  algorithm  of  the  present  paper  in  some  detail.  We  start  with  an 
informal  description,  follow  with  a  more  detailed  one,  and  finish  with  a  complexity  analysis. 

4.1  Informal  Description  of  the  Algorithm 

Below,  we  describe  an  FFT-based  lO^'-order  iterative  algorithm  for  the  solution  of  the 
Lippmann-Schwinger  equation 

00 x)-k 2  [  Gk(x,y)V(y)ip(y)dy  =  k2  [  Gk(x,y)V(y)  <p(y)  dy  (119) 

J  D  J  D 

in  two  dimensions,  where  D  =  [—a,  a]  x  [—a,  a],  Gk  is  the  Green’s  function  for  the  Helmholtz 
equation  in  two  dimension,  i.e.,  Gk(x,y)  =  —  |  •  HQ1\k\\x  —  y||),  and  V(x)  denotes  the 
potential  at  a  point  x.  Here,  ^(x)  and  (j>(x)  are  the  scattered  and  the  incident  fields  at  a 
point  x,  respectively. 

As  discussed  in  Remark  3.9,  once  the  scattered  held  ^  in  the  domain  D  is  known,  the 
scattered  held  -0  outside  D  can  be  calculated  via  the  standard  trapezoidal  rule  applied  to 
(119).  Therefore,  we  focus  on  obtaining  the  solution  of  (119)  for  x  G  D.  Obviously,  (119) 
can  be  written  as  the  linear  system 


(1  -  A)ip  =  A0,  (120) 

where  0  is  the  unknown  scattered  held  in  R,  0  is  the  given  incident  held  in  R,  1  is  the 
identity  operator,  and  A  is  the  integral  operator  in  (119).  As  discussed  in  Section  3,  we 
use  (114)  to  approximate  the  integral  operator  A  on  the  functions  ip,  <p.  With  the  help  of 
FFT  (see  Remark  3.8),  we  apply  the  discretized  version  of  A  rapidly  to  arbitrary  vectors,  and 
solve  the  linear  system  (120)  iteratively.  We  use  one  of  the  most  popular  iterative  solvers, 
Bl-CGSTAB  (the  stabilized  bi-conjugate  gradient  method)  (see  [10],  [12]). 

4.2  Detailed  Description  of  the  Algorithm 

Comment  [Choose  principal  parameters.] 
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Set  the  size  of  the  scattering  structure  to  [—a,  a]  x  [—a,  a]. 

Set  the  initial  position  of  a  point  source  to  (x0,y0)  to  generate  the  incident  field. 

Choose  precision  e  to  be  achieved  for  the  iterative  solver. 

Choose  an  integer  n;  set  h  —  set  the  number  of  nodes  discretizing  a  side  of  the  square 
to  2 n  +  1,  so  that  the  total  number  of  nodes  in  the  discretization  is  N  —  (2 n  +  l)2. 

Choose  the  wave  number  k  for  the  incident  and  the  scattered  fields. 

Construct  a  two-dimensional  sequence  =  —n,...,n  via  the  formula  Vl3  = 

V(ih,jh). 


Stage  1 

Comment  [Construct  the  values  of  the  the  Green’s  function.] 

For  the  user-specified  h  and  k,  calculate  the  correction  coefficients  D0 ,  D i,  D2,  D3,  D 4,  D5 
in  (71)~(76)  via  interpolation  (see  Remark  3.4). 

Construct  the  two-dimensional  sequence  H  via  the  formulae  (110)  -  (112)  on  the  square 
[—2a,  2a]  x  [—2a,  2a],  and  calculate  its  Fourier  transform  using  the  two-dimensional  FFT. 

Stage  2 

Comment  [Construct  the  right  hand  side  of  the  linear  system  (120).] 

For  a  point  source  ( Xo,yo ),  construct  a  two-dimensional  sequence  $  =  {$ij},i,j  = 
—n, ...,  n  for  the  discretized  incident  field  on  the  domain  [—a,  a]  x  [—a,  a]  via  the  formula  (113). 
Construct  the  two-dimensional  sequence  /  =  •  Vl3},  i,j  =  —n, 

As  in  Remark  3.8,  use  the  two-dimensional  FFT  to  calculate  the  Toeplitz  convolution  of 
the  sequences  H  and  /. 


Stage  3 

Comment  [Solve  the  linear  system  using  iterative  solvers.] 

Use  the  iterative  solver  BI-CGSTAB  to  solve  the  linear  system  (1  —  A)  =  A  (f>  to  the 
predetermined  precision  e.  The  multiplication  Aij)  is  done  via  the  combination  of  FFT  and 
the  Toeplitz  convolution  of  the  two-dimensional  sequences  H  and  g,  where  g  =  •  V%3 } , 

i,j  =  —n,  ■  ■  ■  ,  n  with  =  ip(ih,jh )  (see  Remark  3.8). 

The  solution  is  the  scattered  field  at  the  N  discretization  points  in  the  square  [—a,  a]  x 
[—a,  a]. 


Stage  4 

Comment  [Calculate  the  scattered  field  at  any  point  in  the  two-dimensional  plane.] 

Use  interpolation  to  obtain  the  scattered  field  at  any  arbitrary  point  in  the  square  [—a,  a]  x 
[—a,  a],  based  on  the  scattered  field  at  the  N  discretization  points.  As  in  Remark  3.9,  apply 
the  trapezoidal  rule  to  (119)  to  obtain  the  scattered  field  at  any  arbitrary  point  outside  the 
square  [—a,  a]  x  [—a,  a]. 
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4.3  Complexity  Analysis 

A  brief  analysis  of  the  complexity  of  the  algorithm  is  given  below. 

In  stage  1,  the  construction  of  the  two-dimensional  sequence  H  costs  0(N ),  where  N  is 
the  total  number  of  the  discretization  points  on  the  square  [—a,  a]  x  [—a,  a],  i.e.,  N  =  (2n+l)2. 
The  two-dimensional  FFT  costs  0(Nlog(N)).  Thus,  the  CPU  time  cost  of  the  stage  1  is  of 
the  order  0(Nlog(N)). 

In  stage  2,  the  construction  of  the  two  dimensional  sequences  <h,  /  costs  O(N),  and  the 
two-dimensional  FFT  costs  0(N  log(N)).  Thus,  the  CPU  time  cost  of  the  stage  2  is  of  order 
0(N\og(N)). 

The  CPU  time  cost  of  the  stage  3  is  of  order  0(Niter  ■  iVlog(iV)),  where  Niter  is  the 
number  of  iterations  required  by  the  iterative  solver  to  get  the  pre-detcrmined  precision  e. 

In  stage  4,  the  CPU  time  cost  of  interpolating  the  held  at  any  point  in  [—a,  a]  x  [—a,  a] 
O(N). 

Summing  up  the  CPU  times  above,  we  obtain  the  time  estimate  for  the  algorithm 

T  =  a(Niter  ■  N  log(iV))  +  (3  ■  N  +  7,  (121) 

where  N  is  the  total  number  of  discretization  points,  Niter  is  the  number  of  iterations  required 
by  the  iterative  solvers  to  reach  the  precision  e,  and  the  coefficients  a,  /3, 7  are  determined 
by  the  computer  system,  implementation,  etc. 

The  storage  requirements  of  the  algorithm  are  determined  by  the  total  number  of  dis¬ 
cretization  points  N  and  the  number  of  iterations  K  performed  before  restarting  the  iterative 
solvers,  and  are  of  the  form 

S  =  0{K-N).  (122) 


5  Numerical  Examples 

The  algorithm  of  Section  4  has  been  implemented  in  FORTRAN  77  in  double  precision. 
In  this  section,  we  illustrate  the  performance  of  the  scheme  as  applied  to  two  scattering 
objects:  a  Gaussian  and  a  crude  model  of  the  human  skull.  The  experiments  were  carried 
out  on  a  2.8  GHz  Pentium  D  desktop  with  2  Gb  of  RAM  and  an  L2  cache  of  1  Mb.  The 
calculations  reported  in  Tables  1  and  3  were  carried  out  with  a  requested  accuracy  of  1CU13; 
the  calculations  reported  in  Tables  2  and  4  were  carried  out  with  a  requested  accuracy  of 
10-9.  We  restarted  the  BI-CGSTAB  every  5  steps. 

Tables  1-4  illustrate  the  numerical  behavior  of  the  scattered  held  at  arbitrary  far-held 
points,  generated  by  the  potential  V  defined  in  (3);  the  incident  held  is  produced  by  a  single 
point  source.  In  Tables  1  and  2,  we  set  the  potential  V(x,y)  =  e-40^  +y  \  Tables  3  and  4 
illustrate  the  numerical  behavior  of  the  scattered  held,  generated  by  a  model  of  the  human 
skull,  as  shown  in  Figures  2  and  3.  The  headings  of  the  Tables  are  as  follows: 
k  is  the  wave  number  defined  in  (2); 

the  computational  grid  is  A  x  A  for  a  total  of  N2  discretization  points; 
the  computational  grid  is  sizeQbj  wavelengths  x  size0bj  wavelengths; 

N\  is  the  number  of  discretization  points  per  wavelength; 
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Figure  2:  The  human  skull  model 
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Figure  3:  The  human  skull  model  viewed  from  the  top 


Table  1:  10th  order  convergence  of  the  algorithm  for  Gaussian  objects 


k 

N 

size0bj 

Nx 

Erei 

N iter 

tcpu 

25 

8A 

6.28 

6.33E-06 

16 

1.2E-01 

25 

8A 

12.6 

6.63E-09 

16 

5.9E-01 

25 

8A 

25.1 

6.04E-12 

16 

2.6E+00 

25 

400 

8A 

50.2 

7.25E-13 

16 

1.1E+01 

25 

800 

8A 

100 

6.32E-13 

16 

5.5E+01 

25 

1600 

8A 

201 

- 

16 

2.4E+02 

25 


Table  2:  Gaussian  objects  with  a  fixed  number  of  points  per  wavelength 


k 

N 

sizc0f)j 

Nx 

Nrel 

Niter 

tcpu 

25 

8A 

6.28 

6.33E-06 

14 

1.1E-01 

50 

16A 

6.28 

3.80E-06 

20 

7.2E-01 

100 

32A 

6.28 

4.44E-06 

33 

5.2E+00 

200 

400 

64A 

6.28 

8.26E-06 

61 

4.4E+01 

400 

800 

128A 

6.28 

1.60E-05 

171 

6.2E+02 

800 

1600 

255A 

6.28 

- 

891 

1.4E+04 

Table  3:  10th  order  convergence  of  the  algorithm  for  the  simulated  human  skull 


k 

N 

size0bj 

Nx 

Erei 

Niter 

tcpu 

25 

8A 

6.28 

1.18E-04 

134 

1.1E+00 

25 

8A 

12.6 

1.91E-07 

133 

4.7E+00 

25 

8A 

25.1 

2.05E-10 

134 

2.1E+01 

25 

400 

8A 

50.2 

4.56E-12 

135 

9.7E+01 

25 

800 

8A 

100 

7.55E-12 

132 

4.7E+02 

25 

1600 

8A 

201 

- 

132 

2.0E+03 

Table  4:  The  simulated 


human  skull  with  a  fixed  number  of  points  per  wavelength 


k 

N 

sizCobj 

Nx 

Erei 

Niter 

tcpu 

25 

8A 

6.28 

1.17E-04 

97 

7.6E-01 

50 

16A 

6.28 

1.55E-05 

165 

5.8E+00 

100 

32A 

6.28 

1.03E-05 

328 

5.2E+01 

200 

64A 

6.28 

1.69E-05 

756 

5.5E+02 

400 

128A 

6.28 

2.21E-05 

3286 

1.2E+04 

800 

255A 

6.28 

- 

13568 

2.1E+05 
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Erel  is  the  average  of  the  relative  errors  of  the  solution  for  the  scattered  held  at  twenty 
chosen  far-held  points; 

Niter  is  the  number  of  iterations  used  by  the  BI-CGSTAB; 

tcpu  is  the  CPU  time  required  in  seconds. 

The  following  observations  can  be  made  from  the  tables  above,  and  from  the  more  detailed 
numerical  tests  performed  by  the  authors. 

1.  For  smooth  scattering  objects,  the  numerical  algorithm  of  Section  4  displays  10t/l- 
order  convergence;  the  CPU  time  required  to  obtain  requested  precision  is  proportional  to 
Niter  •  N'2  log  N,  where  N2  is  the  total  number  of  discretization  points,  Niter  is  determined  by 
the  requested  precision,  the  number  of  iterations  before  restarting  the  iterative  solver,  the 
size  and  the  structure  of  the  scattering  objects. 

2.  For  sufficiently  smooth  scatterers,  the  relative  precision  of  the  solution  is  determined  by 
the  number  of  discretization  points  per  wavelength.  For  example,  to  obtain  5-digit  precision, 
we  need  roughly  6.5  points  per  wavelength.  Thus,  with  our  constraint  of  2  GB  RAM,  five 
digits  can  be  obtained  for  scattering  objects  as  large  as  300  wavelengths  x  300  wavelengths. 

3.  The  number  of  iterations  increases  dramatically  as  the  size  of  the  scattering  object 
increases,  as  shown  in  Tables  2,  4. 

6  Conclusions 

In  this  paper,  we  construct  an  iterative  algorithm  for  the  solution  of  two-dimensional  forward 
scattering  problems.  The  scheme  is  based  on  the  combination  of  high-order  quadrature 
formulae,  rapid  numerical  application  of  the  integral  operator  in  the  Lippmann-Schwinger 
equation,  and  the  stabilized  bi-conjugate  gradient  method  (BI-CGSTAB).  As  illustrated  via 
several  numerical  examples,  the  scheme  is  (2 p  +  4)th  ( p  =  0, 1,  2,  3, ...)  order  convergent;  the 
computational  complexity  of  the  algorithm  is  0(Niter  ■  N2  log  At),  where  Niter  is  the  number 
of  iterations  used  by  the  iterative  solver,  and  N2  is  the  total  number  of  discretization  points. 

The  approach  we  use  for  the  design  of  high  order  center-corrected  quadrature  formulae 
introduced  in  this  paper  is  not  limited  to  functions  of  the  form  (13)  in  two  dimensions;  it 
is  also  applicable  to  functions  of  the  form  (14)  in  three  dimensions,  as  well  as  many  similar 
situations.  Furthermore,  the  method  does  not  require  access  to  each  of  the  functions  P,  Q 
separately  in  (13)  and  (14);  it  only  requires  the  evaluation  of  the  whole  kernel  K  given  a  pair 
of  points  ( x ,  y).  Quadrature  formulae  of  order  higher  than  10  can  also  be  constructed,  though 
the  derivations  become  more  tedious.  Finally,  the  scheme  is  easily  extended  to  rectangular 
regions  of  the  form  [—a,  a]  x  [—5,5],  even  though  this  paper  only  discusses  on  the  square 
region  [—a,  a]  x  [—a,  a]. 

Acknowledgments.  The  authors  would  like  to  thank  Mark  Tygert  for  helpful  discus¬ 
sions. 
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Appendix  A 


Lemmas  6.1  -  6.5  below  provide  analytical  formulae  for 

(x*“  V-1  H0(ky/ x 2  +  y 2))  dx  dy-  ^  ((p/i)*_1(g/i)J_1  H0{k^/ (ph)2  +  (gh)2))-h2, 

(P,q)  7do,o) 

(123) 

with  (i,  j)  =  {(3, 1),  (5, 1),  (3,  3),  (7, 1),  (5,  3)}.  The  proofs  are  straightforward  and  tedious, 
and  use  the  help  of  Mathematica. 

Lemma  6.1  For  any  k  G  C+  and  h  >  0, 


r»oo  /»oo 


D 1  = 


'  — oo  — oo 


dA 


7/q1^  ( kr)x 2  dxdp  —  E  H(o\k\f  ( ph )2  +  (gh)2)  •  (ph)2  ■  h2 

0,0) 

h4  eiaih(eiai/l  +  l)(eiQ2/l  +  1) 

~~  ~2~ 


7T  i-oo  V&2  -  A2  V (7«i)3(i«2) 
4  f”  dA  2 


ia\h  _  ^3^ia2h  _  ^ 


7T  7-oo  Vk2  -  A2  afa2(l  +  Xi)3(l  +  x2) 


+  3^i  +  +  3  y\  +  -pi(iaih)  +  3pip2  +  -yi{ia2h)  +  -y2{iaih) 


( (ogh)2 

V  12 

+x'l  +  3x1X2  +  x^x2  -  -ia2hy2  +  -a\K2{x\  +  2xi) 

3  1 

+-a1a2h2(x1x2  +  Xi  +  z2)  +  +  Zi)2(l  +  z2)  | , 


(124) 


(125) 


where  r,ai,a2,xi,x2,  yi,y2,  zi,z2  are  defined  by  (88),  (89),  and  (93)-(95). 


Lemma  6.2  For  any  k  G  C+  and  h  >  0, 


f*  oo  /*oo 


where 


Do  = 


'  — oo  — OO 


H^\kr)  x4  dx  dy  —  ^  H^\k\J (; ph )2  +  (g/i)2)  •  (p^)4  •  h2 

(p,q)?(  o,o) 


dA 


7T  7_oo  \/fc2  -  A2 

24  ^6  giaih  _j_  4^2*01/1  -)_  He^*ai^  +  e4iaih  e*“2fe  +  1 

2" 


(iai)5(m2) 

4  f00  dX 


_  4^5 


gia2  h  4 


7r  V&2  -  A2  afa2(l  +xi)5(l  +  x2) 

3  20 

2  (ia2/i)2  +  (ia2h)3  +  5(iai/i)(ia2/i)2  +  —  (ia2h)4  +  —  (iaih)2  (ia2h)2 

1 U  O 

5  1  45 

+-(iaih)(ia2h)3  +  -{iaih)4(ia2h )  +  4a;1(ia:1/i)4  +  —  x^iiaih)3  (ia2h) 

+2xi(iaih)4  {ia2h)  +  A3x\{ia\h)3  +  25o;2(ia;i/i)2(ia:2/z)  +  Qx\(ioi\h)4 
45 

+— x\(iaih)3  (ia2h)  +  3a;2(iai/i)4(ia:2/i)  +  lbx\(iaih)3  +  Ax\{ia\h )4 

15  1 

+— x\{ia.ih)3  (ia2h)  +  2xf  (iaih)4(ia2h)  +  x4(iaih)4  +  -x4(ia:i/i)4(ia;2/i) 

15  1 

+— x2(ia\h)3  (ia2h)  +  - x2(iaih)4(ia2h )  +  50xia;2(iQ;i/i)2(ia2/i) 

45 

+— XiX2{iaih)3  {ia2h)  +  2xix2(icxih)4{ia2h)  +  25x2a;2(ia:i/z)2(ia;2/i) 

45  15 

+— x\x2(iaih)3  (ia2h)  +  ?>x\x2(iaih)4{ia2h)  +  —  x\x2{ia\h)3  {ia2h) 

-\-2x\x2(ioiih)4(ioi2h)  +  -£4a;2(ia:i/i)4(m2/i)  +  12w2(ia2h)  +  30z2(ia:i/i)(ia2/i) 

+40y2(ia;i/i)2(ia:2/i)  +  I5yi(iaih)  (ia2h)2  +  30yiy-2(iaih)  (ia2h)  —  120«i  —  24m2 
— 180wi(iai/i)  —  60w2(iaih )  —  60wi(io;2/z)  —  Vd§Zi(ioiih)2  —  90zi(ia!i/i)(*a!2/i) 
— 80x:2  (icti/?)2  +  35yi(ia\h)3  —  AOyfiiaih)2  (ia2h)  —  I30yl(iaih)2 

—  180yl(ia?ih)(ia2h)  —  2AQy\{iaih)  —  120  y\{ia2h)  —  120y4  —  30y2(iaih)3 

—  180yiy2(iaih)2  —  360yly2(iaih,)  —  240 y\y2  —  2Qzi(ia2h)2  —  240^2  —  120ziz2 
-360y2(iai/z)  -  240 y\  -  120yf(ia2h)  -  240t/i7/2  (*04/1)  -  240 y\y2 

— 2\.x\  —  120x4a;2  —  2Ax\x2  \ , 


Wi  =  z  i 


U\  =  W\  — 


(■ ia\h )3 
24 

(iai/i)4 

120 


( iaih)n 

§(^’W!  = 


^2 


n=4 

oo 


>.  7 — rTAT’^2  =  W2  - 


(■ ioL2h )3 
24 

0 ioi2h )4 
120 


(7a2/i)n 
^  (n+T)!’ 

(ia2/i)n 


E 

n=4 

oo 

E 


n=5  v  7  n= 5 

and  r,oii,ot2,Xi,x2,yi,y2,Zi,z2  are  defined  by  (88),  (89)  and  (93)-(95). 


ri(n+  1)!’ 


(126) 

(127) 

(128) 
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Lemma  6.3  For  any  k  G  C+  and  h  >  0, 


D,  = 


H^\kr)  x2  y2  dx  dy  —  ^  Hj)1\ky/ ( ph )2  +  ( qh )2)  •  ( ph)2(qh )2  •  h 2 

(p,  9)^(0,  o) 


0iaih Sgiaih  _j_  -Q  giaift^giaift 
^giai/i  _  1)3  ^giaih  _  1)3 


7T  ^/A;2  -  A2  V(*«l)3(*Q;2)3  (eiaih  _  1)3  (eiaife  _  l)3  / 

-  1  Z00  d\  _ -4 _ 

7T ./  oc  \/A;2  -  A2  a?a^(l +a;i)3(l +  x2)3 

'  (2I) ' ^ai/l^  +  2io(-2^)4  +  ^i(iai^)2  +  |-2(*a2^)2 

9  15  15  3 

+-(iaih)(ia2h)(zi  +  z2)  +  —  z1(ia2h)2  +  — £2(*ai/*)2  -  -yi(ia:i/i)2(ia:2/*) 

3  9  9 

-- y2(iaih)(ia2h)2  +  3m  1  +  3u2  +  -w2(iaih)  +  -Wi(ia2h) 

3  3  9 

+-Wi(io;i/i)  +  - w2(ia2h )  -  —  (zcki /i) (ia;2/i.) (^1  +  z2)  +  3^2  +  3 z\ 

1111 

+9^1^2  -  -yii^iaihY  -  - y2(ia2h )3  -  -yf(ia i/i)2  -  -1 y\{ia2h )2 

9  119 

+-(*Q!i/i)(iQ!2/i)(3j/ij/2  +  -yi(iot2h)  +  -y2(*a!i/i))  +  -yi(iai/*)2(ia2/i) 

4  008 

3  9  9.3 

+  ^2/i(4Q!2/i)3  +  -y2(iai/z)(io!2/i)  +  - -y\[ia2h )2  +  - y\(ia2h ) 
o  4  4  z 

3  9  9  9 

+  -y2(iai/z)3  +  -y2{ioi\h)(ia2h)2  +  -y^fiiaih)2  +  -yiy2(ia2h)2 
00  4  4 

9  9  9 

+  2^2(*ai^)  +  9  y\y2(ia.2K)  +  3  y3y2  +  -yl(iaihf  +  -j/|(ia:i/i)(iQ!2/i) 

9  3 

+9yiy|(ia;i/i)  +  ^ywl{ioi2h)  +  9y2y2  +  -y2{iot\h)  +  3j/iy| 

— -(iai/z)2(io;2/z)2(2a;i  +  x2  +  2a;2  +  4xia;2  +  2x\x2  +  £2 

3  3  9 

+2x^1  +  x2^)  +  -y2(ia  1/1)  +  y3  +  -yl(ia2h)  +  y2  +  - y\[ia2h ) 

9 

+9yiy2(ia;i/i  +  *0(2/1)  +  9y2y2  +  -y2(ioi\h)  +  9yiy2 
3 

--(*a:i/*)2(ia:2/*)(a;2  +  2;cia;2  +  xfx2) 

■~-(ia1h)(ia2h)2(x2  +  2x1x2  +  x\x\)  +  3x3a;2  +  3x2a;2  +  x\x\\ , 


(129) 


(130) 


where  r,  oti,  a2,xi,x2,yi,y2,zi,z2,wi,w2,ui,u2  are  defined  by  (88),  (89),  (93)-(95),  (127) 
and  (128). 
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Lemma  6.4  For  any  k  G  C+  and  h  >  0, 


'•oo  poo 


Dd  = 


/  — oo  J  — oo 
r*QO 


Hq1'1  ( kr )  x 6  dx  dy  —  ^  (ky/ ( ph )2  +  (g/i)2)  •  (p/i)6  •  /?2 

(P, 9)^(0, 0) 


dA 


7T  7_oo  V&2  -  A2 

/  720  /z8  eiai/l  +  57e2iQl/l  +  302e3iQl/l  +  302e4iaih  +  57e5iQl/l  +  e6iaih  eia2h  +  1 ' 

\ala2  2  ( eiaih  -  l)7  eia^h  - 

4  f°°  dX  1 

vr  L  Vk2  -  A2  aja2(l  +  £i)7(1  +  x2)  ^  + 


where 


A  =  60  (a2h)2  —  30  (ia2h)3  —  210(iaih)(ia2h)2  —  9(a2h)4  —  105(aih)(a2h )3 
-385(ai/z)2(a2/i)2  +  37800j/?(iai/i)  +  25200 y\  +  7560j/?(ia2/i) 

+15120y1y2(ia1h)  +  15120  yfy2  +  15120x5  +  5040x5  +  720x[  +  25200^2 
+15120x5x2  +  5040x®x2  +  720xjx2  +  5040«i  +  720m2  +  I2600wi(ioiih ) 

—360 w2(ia2h)  +  2520w2(io;i/i)  +  2520w  i(ia2h)  —  17220(a;i/i)2£i 
+6300^1  (iaih)(ia2h)  +  I260z2(aih)(a2h)  —  4620(o;i  h)2z2  +  15120^2 
+2940yi(io;i/z)3  —  Q999y\(a\h)2  (ia2h)  —  34440y2(«i/i)2  —  18900yl(aih)(a2h) 
+50400y3(mi/i)  +  12600?/3(ia2/i)  +  25200^  +  3150y2(iai/i)3 
— 18900yi|/2(cn^)2  +  37800y2y2(io;i/i)  +  25200 y\y2  —  8A9zi(a2h)2  +  5040xi£2 
+2310y2(o;i/i)2  (ia2h)  +  630i/i(iai/i)(a:2/i)2  +  1260t/i7/2(q;i/i)(q;2/i) 

63  1 

— 63(ia.ih)5  —  301  (iaih)4  (ia2h)  +  (cti h)6  +  —  (aih)5  (a2h)  +  -(aiKf  (ia2h) 

— 2408xi(o;i/i)4  —  3150xi(ai/i)3(a2/i)  —  315xi(ia:i/i)5  —  1204xi(iai/i)4(ia:2/i) 

315 

+6xi(ai h)6  H — —xi(aih)5 (a2h)  +  3xi(ai h)6(ia2h)  —  6300x5 (i«i^)3 

+1680x2(ai/z)2(ia2/z)  —  3612x2(q;i/i)4  —  3150x2(ai/z)3(o;2/z)  —  630x2(i«i/z)5 

15 

— 1806x2(iai/i)4(ia2/i)  +  15x2(q;i  h)6  +  315x2(q;i/i)5  (a2h)  +  — x2(cti  h)6(ia2h) 

— 2100x5(iai^)3  —  2408x5  (ax/?)4  —  1050x5(ai/z)3(a2/z)  —  630x5  (iaih)5 
— 1204x5  (iaih) 4  (ia2h)  +  20x5(ai/i)6  +  315x5(ai/i)5(a2/i)  +  10x5(ai/i)6(ia2/i) 

— 602x5(o;i^)4  —  315x4(iaih)5  —  301x4  (iaih)4  (ia2h)  +  15x5(ai/z)6 
315  15 

H — —x4(ai  h)5(a2h)  +  —  x\(aih)Q  (ia2h)  —  §3x\(iaih)b  +  6x5(ai/z)6 

63  1 

+— xl(aih)5  (a2h)  +  3x\(aih)& (ia2h)  +  x\(aih)G  +  -Xi(otih)6  (ia2h) ,  (132) 
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B  =  —  1050a;2(a;i/})3(a2^)  —  301x  2(aih)4  (ia2h)  +  —  x2(aih)b  (a2h) 

-\--x2{aiKf,{ia2h )  +  3360a;i£2(a:i/i)2(ia:2^)  —  31§Qx\x2(oi\h)3  (a2h) 

315 

—  l204xiX2(aih)4(ia2h)  H — —  x1x2(a.ih)b  (a2h)  +  3x1x2(a1h)6(ia2h) 

+  1680xlx2(aih)2  (ia2h)  —  3150xlx2(aih)3  (a2h)  —  1806x2a;2(cti^)4(ia;2^) 

15 

+313x\x2{a\Ky  {a2h)  +  —  x\x2(oi\h)& {ia2h)  —  lD3Dx\x2{a\h)3  (a2h) 

—  120Axlx2(aih)4(ia2h)  +  315xlx2(aiti)5  (a2h)  +  I0x3x2(aih)e(ia2h) 

315  15 

— 301x)x2(aih)4  (ia2h)  3 — —x)x2(aih)5  (a2h)  +  —  x4x2(aih)6  (ia2h) 

63  1 

+— x\x2{oiih)b  (a2h)  +  3x\x2(oiiKf  (ia2h)  +  -x\x2(oiiKf  (ia2h) ,  (133) 

and  where  r,a\,a2,xi,x2,yi,y2,  z\,  z2,wi,w2,ui,u2  are  defined  by  (88),  (89),  (93)-(95), 
(127)  and  (128). 


Lemma  6.5  For  any  k  G  C+  and  h  >  0 
D5  = 


r»oo  /»oo 


'  — oo  J  — oo 


d\ 


7T  Vk2  -  A2 
4  dX 


H^fikr)  x4  y2  dx  dy  —  ^  H^\k\/  ( ph )2  +  ( qh )2)  •  (ph)4(qh)2  ■  h2 

(P, 9)^(0, 0) 

4g  ^  giai/i  _J_  j^g2*alfe  _|_  -j^g3iai/i  g4*“lfe  .  ^g*a2/i  _|_  ^ 


afa3 


7T  7-00  a/A:2  -  A2  afa^(l +Xi)5(l +  a;2)3 


•(C'  +  D) 


_  ]^3 


(134) 


where 

l 

C  =  ~(a2h)4  +  240wi  +  360wi(io;i/z)  +  144m2  +  360w2(mi/4) 

5 

+360wi(ia2^)  +  72w2(ia2h)  —  320z\(aih)2  —  3AQzi(aih)(a2h) 

—A8Qz2(oiih)2  —  300^i  {a2h)2  —  180z2(a?ih)(a2h)  —  36z2(a2h)2  +  48CL2 
+720^1^2  +  144^2  —  70yi(ia?ih)3  —  2A0yi(aih)2  (ia2h)  —  210yi(iaih)(a2h)2 
+30yi(ia2h)3  —  260y2(ai/i)2  —  1080yl(aih)(a2h)  —  360 y((oi2h)2 
+A8Qy\[iaih)  +  72tiy\{ia2h)  +  240y4  +  18Qy2{ia\h)3  +  120y2(iaih)2  (ia2h) 
+30?/2(*ai/z)(Q;2/4)2  -  1080?/i|/2(ai^)2  -  1260yiy2(ai^)(a2^) 

-180yi|/2(a2/i)2  +  2imy2ly2(ioilh)  +  U40yly2(ia2h)  +  144%il/2  -  360y2(a1h)2 
—  l8Qyl(aih)(a2h)  +  144:0yiy2(iaih)  +  360yiy2(ia2h)  +  1440 y\y\  +  l2Qy2{iaih) 
+240?/i|/2  +  72Qy\{iaih)  +  480y3  +  720y((yia2h)  +  1440?/i|/2(*Q:i^) 

+1440t/27/2  +  720y1y2(ia2h)  +  360y2(ia1h)  +  720y1yl  +  72y\{ia2h)  +  A8y32  +  48x'i 
+720a;4a:2  +  1AAx\x2  +  1440;r3a;2  +  720x)x2  +  1AAx^x2  +  480a;2a;2  +  480o;32;2,  (135) 
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D  =  240x^2  +  48x^2  —  24y2(ia2h)3  +  24yl(a2h )2  —  3{ia\h)4{iot2h ) 

— 15(iaih)3  (ia2h)2  +  (aih)4(a2h)2  —  8xi(q;i  h)4  —  I33x  fiaih)3  (a2h) 

—  IWxfiaih)2  (a2h)2  —  I2xfia.ih)4(ia2h)  —  A5xi(iaih)3  (ia2h)2  +  Axi(aih)4(a2h)2 
— QQx\{ioiih)3  —  lbQx\{iaih)2  (ia2h)  —  I2xl(a±  h)4  —  135xl(a?ih)3  (a2h) 

— 50xl(a?ih)2  (a2h)2  —  18xl(aih)4  (ia2h)  —  A5x\{ia.ih)3  (ia2h)2  +  6xl(aih)4  (a2h)2 
—30x3  (iaqh)3  —  8x\(a.\h)4  —  A3x\(aih)3  (oi2h)  —  I2x\(a.\h)4{ia2h) 

— I3x\(iaih)3  {ia2h)2  +  Ax\(aih)4(a2h)2  —  2x\(a1  h)4  —  3x\(a.ih)4(ia2h ) 
+x4(ai/z)4(o!2^)2  —  4:5x2(aih)3(a2h)  —  lOOxfiafii)2  (a2h)2  —  3x2(aih)4  (ia2h) 

— 30x2(iaih)3  (ia2h)2  +  2a;2(ctih)4(a2h)2  +  300xia;2(ai/i)2(ia;2^) 
+I20xix2(ia.ih)(a2h)2  —  13hx\X2{oi\h)3  (a2h)  —  200xix2(aih)2  (a2h)2 
—12x\x2(a.ih)4(ia2h)  —  99xix2(ia\h)3  (ia2h)2  +  8x1x2(a\h)4(a2h)2 
+I39x2lx2(aih)2  (ia2h)  —  135a;2a;2(cn^)3(a2^)  —  100£2£2(a:i/})2(a;2^)2 
—I8x\x2(aih)4{ict2h)  —  99x\x2{iatih)3  (ia2h)2  +  ^^^(ctih^a^h)2 
— 45x3a;2(ai^)3(a2^)  —  I2x\x2(aih)4{ia2h)  —  39x\x2[iaih)3  (ia2h)2 
+8x\x2(o>ih)4(a2h)2  —  3x4x2(aih)4(ia2h)  +  2x4x2(o;i/i)4(q;2^)2 
+60xl(iaih)(a2h)2  —  50xl(aih)2  (a2h)2  —  15xl(iaih)3  (ia2h)2  +  x22(a\h)4  (a2h)2 
+Q9xix\(iaih)(a2h)2  —  lOOxia^aqh)2  (a2h)2  —  A3xix\{iaih)3  (ia2h)2 
+Axixl(oiih)4(a2h )2  —  39x\x\{aih)2  (a2h)2  —  kf>x\x\{ia.ih)3  (ia2h)2 
+Qxlxl(aih)4  (a2h)2  —  l§x\x^(iaih)3  (ia2h)2  +  Ax\x22(a\h)4  (a2h)2 
-\-x\x^(aih)4(a2h)2 ,  (136) 

where  r,  an,  a2,  x\,  x2l  y±,  y2,  z\,  z2,  w\,  w2,  u\,  u2  are  defined  by  (88),  (89),  (93)-(95),  (127) 

and  (128). 


Appendix  B 


Here,  we  present  the  center-corrected  quadrature  formulae  of  orders  4,  6  and  8  for  the  integral 


4>(x,  y)  ■  Hfik\Jx2  +  y2)  dxdy. 


(137) 


Lemma  6.6  (The  4th-Order  Center-corrected  Quadrature  Formula).  Suppose  that 
n  >  1  is  an  integer,  and  a,  h  are  two  positive  real  numbers  such  that  h  =  a/n.  Suppose 
further  that  0  :  R2  — >  C  is  a  function  such  that  <fi  G  c2  (R  x  R) ,  and  that  <p  is  zero  outside 
the  square  [—a,  a]  x  [—a,  a].  Then,  for  any  k  G  C+, 


4>(x,  y)  ■  Hfiksjx2  +  y2)  dx  dy 


UTh  (f(x,  y)  ■  H0(k^/x2  +  y2)')  +  0{h4). 


(138) 


In  (138), 

UTh[fi(x,y)-H0(ky/^+7))  =rh^(x,y)-H0(k^T72))  +  ^  ^(ph,  qh),  (139) 

p,qeS 
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where 

S  =  {p,  q  E  Z  :  p  =  0  and  q  =  0},  (140) 

T'h  (^(x,  y)  ■  Ho(k^Tf))  =  fab’  *h)  '  Ho(ky/(ph)>  +  ( qh )2))  •  h2,  (141) 

(P, 9)^(0, 0) 

and 

roo  =  D0.  (142) 


Lemma  6.7  (  The  6^-Order  Center-corrected  Quadrature  Formula).  Suppose  that 
n  >  1  is  an  integer,  and  a,  h  are  two  positive  real  numbers  such  that  h  =  a/n.  Suppose 
further  that  <f  :  R2  C  is  a  function  such  that  (j)  E  c4(R  x  R),  and  that  <f  is  zero  outside 
the  square  [—a,  a]  x  [—a,  a].  Then,  for  any  k  E  C+, 


'  —a  J  —a 


(j){x,  y)  ■  H0(ky/x 2  +  y2)  dxdy  =  UTh  y)  ■  H0(k\/x2  +  y2)  j  +  0(h6).  (143) 


In  (143), 

UTh{f)(x,y )  ■  H0(k\/x2  +  y2))  =  T,h(^(x,y)  ■  H0(ky/ x2  +  y2)^  +  ^  r^plpqh),  (144) 

p,qeS 


where 

S  =  {p,  q  E  Z  :  \p  +  q |  <  1  and  \p  —  q\  <  1}, 

T'h  (f>{x,  y)  •  H0(k^T72))  =  fa11’  <*h)  '  Ho(W(ph)2  +  (qh)2))  ■  h2 , 

(p,  q)^(0, 0) 


and 


roo  ~  2 


T±10  —  r0±l 


£i 

h 2’ 

m 

2  /i2  ' 


(145) 

(146) 


(147) 

(148) 


Lemma  6.8  (  The  8*h-Order  Center-corrected  Quadrature  Formula).  Suppose  that 
n  >  1  is  an  integer,  and  a,  h  are  two  positive  real  numbers  such  that  h  =  a/n.  Suppose 
further  that  <j>  :  R2  C  is  a  function  such  that  f>  E  c6(R  x  R),  and  that  <f  is  zero  outside 
the  square  [—a,  a]  x  [—a,  a].  Then,  for  any  k  E  C+, 


4>(x,  y)  ■  Hq  ( k \Jx?  +  y2)  dx  dy 


UTh  y)  ■  H0(k  V^T^))  +  0(h8). 


(149) 


In  (149), 

UTh(j>(x,y)  ■  Ho(k^x2  +  y2)^  =  T'h(j>(x,y)  ■  H0(k^/x2  +  y2)^  +  ^  Tpq(j>(ph,qh),  (150) 

p,qeS 


where 

S  =  {p,  q  E  Z  :  \p  +  q\  <  2  and  \p  —  q\  <  2},  (151) 
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T'h(f>{x,y)  ■  i70(fcvWs/2))  =  E  (<K^^)  '  H0(W(ph)2  +  W))  '  ^2,  (152) 


(p,«)  tRo,o) 

h  n  5  D1  1 D2  D3 

00  0  2  h2  +  2  h4  +  h4  ’ 

(153) 

h _ h  _  2  D1  1  D2  1  D, 

±10  o±i  3  fpi  q  h\  2  h4  ’ 

(154) 

h  _  h  _  1  Di  1  D2 

±20  0±2  24  h2  +  24  IP  ’ 

(155) 

Th  _  1  D3 

±1±1  4  IP  ' 

(156) 
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